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Abstract 

The software tool GRworkbench is an ongoing project in visual, 
numerical General Relativity at The Australian National University. 
Recently, GRworkbench has been significantly extended to facilitate 
numerical experimentation in analytically-defined space-times. The 
numerical differential geometric engine has been rewritten using func- 
tional programming techniques, enabling objects which are normally 
defined as functions in the formalism of differential geometry and Gen- 
eral Relativity to be directly represented as function variables in the 
C++ code of GRworkbench. The new functional differential geometric 
engine allows for more accurate and efficient visualisation of objects in 
space-times and makes new, efficient computational techniques avail- 
able. Motivated by the desire to investigate a recent scientific claim 
using GRworkbench, new tools for numerical experimentation have 
been implemented, allowing for the simulation of complex physical 
situations. 

1 Introduction 

Physically important exact solutions of the Einstein field equation 
are often difficult to work with algebraically because of their com- 
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plexity. It is usually necessary to make simplifying assumptions and 
approximations if analytic results are desired. The goal of the ongoing 
GRworkbench project at The Australian National University is to cre- 
ate a visual software tool for numerical General Relativity. Working 
with S. M. Scott and B. J. K. Evans, A. C. Searle implemented a new 
version of GRworkbench in 1999. It featured an imbedded platform- 
independent GUI (Graphical User Interface), a novel numerical differ- 
ential geometric engine, and a flexible visualisation system, and it was 
easy to extend with additional space-time definitions. P 

In the highly general visualisation system of GRworkbench, space- 
times are visualised by transforming the 4 coordinates of any space- 
time chart under arbitrary distortions down to a 3-dimensional vi- 
sualisation space, which is rendered on the screen from an arbitrary 
viewpoint using the OpenGL graphics library. 

The differential geometric engine of GRworkbench allows for ab- 
stract objects, such as points of a manifold and tangent vectors, to 
be associated with multiple numerical representations, corresponding 
to the 'image' of the abstract object in different coordinate charts. 
As a part of the definition of each space-time, GRworkbench is sup- 
plied with the maps between the various coordinate systems. Numer- 
ical operations, such as integration of the geodesic equation ('geodesic 
tracing'), are performed in the coordinates of a single chart, until a 
chart boundary or other obstacle is encountered, at which point the 
algorithms transform the data into another coordinate system and 
attempt to resume the computation there. 

Two pieces of information define a space-time in GRworkbench: (1) 
a set of coordinate systems and the maps between them; and (2) func- 
tions giving the components of the metric tensor on each chart. For 
numerical operations, such as geodesic tracing, which involve deriva- 
tives of the metric components, numerical methods are employed to 
compute the derivatives. 

1.1 New developments 

The developments described in this article were motivated in part by 
a desire to employ GRworkbench in an analysis of a recent scientific 
claim. The specific goal was to extend GRworkbench to facilitate the 
simulation of interesting and potentially complex physical situations 
in numerical experiments. 

In Section|2lwe describe how the numerical and differential geomet- 
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ric aspects of GRworkbench have been rewritten and extended using 
functional programming techniques. Functional programming enables 
a more direct representation of those aspects of numerical analysis and 
differential geometry that are naturally defined in terms of functions or 
operations on functions. The functional representation of differential 
geometric aspects in GRworkbench has enabled new computational 
techniques to be developed and has allowed for more accurate and 
efficient visualisation methods to be employed. 

In SectionEJwe describe new tools for modelling physical systems, 
which were implemented to facilitate the investigation of the motivat- 
ing scientific claim. A summary of the original claim and our investi- 
gation of it are described elsewhere. [2j 



1.2 Definitions and notation 

A chart is an open subset C C M", representing a coordinate system 
on a subset Aic C of the n-dimensional space-time manifold A4 
(n = 4). We denote by (pc ■ M.c — > C the one-to-one and onto function 
which maps points in into the chart C. The tangent space of a 
point p € is denoted by Tp. The components of the metric tensor 
g: Tp xTp ^M. are denoted by gij . 



2 Functional differential geometry 

Consider the world- line of a freely- falling particle. Mathematically, 
it is a function /: R — > In GRworkbench , a point p G is 

represented numerically by its coordinates = (t)c{p) ^ I^" on one 
or more charts C . On a chart C the coordinates of the world- line are 
functions of the world-line parameter t. For the freely- falling 

particle, the are the solutions of the geodesic equation^ 

(Px^ dx" dx^ _ 

Previously in GRworkbench, a function such as / was represented 
by its images f{tj) at a finite number of values, tj, j = 1,...,N, 
where typically the tj were evenly spaced: tj = jAt. Using an ode 
integrator from a standard library of numerical algorithms, GRwork- 
bench numerically integrated the geodesic equation to determine the 



^Throughout this article we assume that geodesies are always affinely parameterised. 
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coordinates x^{tj) (on some chart) of the points f{tj). The number 
A'^ of i-values and the separation Ai between them were specified by 
the user according to their requirements; for example, so as to enable 
the world-line to be smoothly displayed in the visualisation system 
of GRworkbench. Such a representation was adequate for the simple 
visualisation and numerical experimentation tasks to which GRwork- 
bench was first applied. 



2.1 Operations on functions 

Consider the operation of parallel transport which, given an initial 
tangent vector f at a point p and a curve / passing through p, uniquely 
determines a tangent vector at each other point on /. On a chart C the 
components of v are determined as functions of the curve parameter 
t by the parallel transport equation 

where the x°'{t) are the coordinates of the curve / on C. 

In order to numerically integrate Equation ^ to any desired pre- 
cision, it is necessary to be able to determine the quantities dx"' /dt for 
any value of t. (It is sufficient for the functions x"(t) to be defined for 
any value of t, from which the dx"'/dt can be obtained via numerical 
differentiation.) This presents no difficulty if the curve / is prescribed 
explicitly, for example, 



x\t) 



t, ifi = 0; 
0, otherwise. 



If, however, the curve is represented only by its coordinates at a finite 
number of values of t, as described above, then there will not, in 
general, be sufficient information to numerically integrate Equation ^ 
to any desired precision. 

Therefore, if we wish to be able to use numerically determined 
curves, such as geodesies obtained by the numerical integration of 
Equation as inputs to algorithms such as the one for numerically 
integrating the parallel transport equation, then we must represent 
those curves in a stronger way than by their coordinates at a finite 
number of points. The cornerstone of functional programming is the 
ability to store operations, such as the numerical integration of the 
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geodesic equation, in program variables, so that, like any other vari- 
able, they can be passed as arguments to functions or created and 
returned as the result of functions. 

2.2 Representing functions directly 

In the functional differential geometric framework of GRworkbench, 
a curve in space-time is represented by an algorithm which, given 
any value of the curve parameter t, computes and returns a point 
p (z A4. For example, the algorithm representing a geodesic / defined 
in terms of initial conditions at t = is one which, given a value tf, 
numerically integrates Equation from t = to t = tf to determine 
the coordinates of f{tf) on some chart, and returns the point which 
has those coordinates. 

In the following sections we describe the advantages of the func- 
tional representation of space-time curves for the two most impor- 
tant utilities of GRworkbench: computation in, and visualisation of, 
analytically-defined space-times. The details of the implementation of 
functional programming techniques in the C+-|- code of GRworkbench 
are described elsewhere. |j3j 

2.3 Application to visualisation 

Space-times are visualised in GRworkbench in a three-dimensional 
space which is rendered onto the two-dimensional computer screen 
from an arbitrary camera position using the OpenGL graphics library. 
The n-dimensional coordinate system of a chart C undergoes an ar- 
bitrary user-specifiable transformation v: down to three 
dimensions for visualisation, the simplest of which is simply the sup- 
pression of all but three of the n coordinates. 

When a space-time curve / : R ^ is visualised on a given co- 
ordinate chart C, GRworkbench is effectively required to render a 
parametric plot of the function 

g-.R^M.^ g = vo(l>cof, (3) 

where, recall, (j)c'- M.c C, M.c C M., is the map associated with 
the chart C C M". The function g is only defined for values of the 
parameter t of f such that f{t) S A4c- In GRworkbench, a plot of 
the function g in consists of straight line segments connecting the 
points g{t) for an increasing sequence of values of t. As more points are 
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Figure 1: A continuous curve is smoothly visualised in GRworkbench by 
ensuring that the angle between any two consecutive straight line segments, 
which comprise the visualisation, is less than a small angle 6. For the curve 
depicted, {n — ZACE\ > 6. After subdivision of the lines AC and CE, 
Itt - ZABC\ < 5 and [tt - ZCDE\ < S. 

added, the plot represents g more accurately. The number of values 
of t at which g (and hence /) must be evaluated depends on how 
accurately we require the curve to be visualised. 

As a concrete example, consider a circular orbit about the origin 
in a spherically-symmetric space-time. In spherical polar coordinates 
{t,r,9,(f)), with curve parameter s, the world- line satisfies 

t = s, r = R, 9 = tt/2, (p = ms, 

where to = d(j)/dt is the constant angular speed and i? is a constant. 
When a segment of this curve is visualised in spherical polar coordi- 
nates, plotting any 3 of the (t, r, 9, 0) coordinates orthogonally (and 
suppressing the remaining coordinate), the result is a straight line. 
Thus, the entire circular orbit can be visualised accurately after eval- 
uating the coordinates of the world-line at only two values of s; a 
straight line joining two points of a linear function is the best possible 
visualisation of that function. 

Suppose, on the other hand, that we wish to visualise a segment 
(from s = a to s = 6) of the same circular orbit in rectangular co- 
ordinates {t,x,y,z), where x = r sin cose/), y = r sin sine/), and 
z = rcosO. The world-line is curved in these coordinates, and thus 
the quality of the visualisation will depend on the number of points 
at which the function is sampled. 

In GRworkbench, the quality of the visualisation of curves is pa- 
rameterised by an angle 6: the angle between any two consecutive 
straight lines in the visualisation of a curve is required to be less than 
6. This requirement is depicted in Figure E 

To satisfy this 'smoothness' requirement while plotting a curved 
function such as the circular orbit under consideration, GRworkbench 
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employs a recursive algorithm. The function ^r: M ^ is evaluated 
at an initial set of points, which may contain as few as three points (the 
endpoints and midpoint of the range of values over which the curve 
is to be visualised). Next, the angle between each successive pair of 
straight line segments is determined and, if the angle is not less than 
6, each line segment is broken into two line segments, each being half 
the length (in the curve parameter) of the original line segment. Each 
such subdivision requires the evaluation of (7 at a new value of the 
curve parameter. 

There are two clear advantages of this method of visualisation 
over the original approach, in which the user was required to specify 
how frequently the function g was to be sampled. Firstly, the user is 
no longer required to know (or guess) in advance how many straight 
line segments will be needed to smoothly visualise a curve. Secondly, 
in regions of higher curvature,^ the function is automatically sam- 
pled more frequently in order to produce a consistently smooth plot, 
and in regions of lower curvature the function is sampled less fre- 
quently, preventing unnecessary computation. This latter advantage 
is an example of so-called 'lazy' function evaluation, and can result 
in a noticeable increase in graphical performance for functions which 
are expensive to evaluate, such as numerically integrated geodesies or 
more complicated objects such as the parallel curves of ^2A\ 

2.4 Application to computation 

The operation of geodesic tracing (numerical integration of Equa- 
tion from initial conditions) can be thought of as a function from 
tangent vectors to space-time curves: 

geodesic: Tp^ (R^ M), 
geodesic('t;) = A, XiR^M, 

A is the geodesic with tangent vector v at p = A(0), (4) 

where, for any sets A and B, we denote by (A ^ B) a set of functions 
from A to B. 

As explained above, an important motivation for adopting func- 
tional programming techniques in GRworkbench was the desire to be 
able to use numerically determined functions, such as geodesic (f), for 

^By 'curvature' here we mean the curvature of the function g: M. ^ of Equation 
not the space-time curvature. 
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some tangent vector v, as operands to other numerical operations, 
such as the parahel transport of a vector along a curve. From the 
point of view of an algorithm for integrating the parallel transport 
Equation ((21), it is not important how the function / is defined: as 
long as / takes a real number as an argument and returns an ele- 
ment of Ai, the algorithm can proceed. Thus, the parallel transport 
operation can itself be thought of as a function: 

paralleLtransport : {R ^ M) xV^{R^V), 
paralleLtransport(/, f ) = h, /i: M — > y, 

h(t) is the parallel transport of f G ^/(o) to f{t) along /, (5) 

where V is the set of all tangent vectors on the space-time mani- 
fold A4. If the curve / is the geodesic world-line of a freely- falling 
non-rotating observer, and v is orthogonal to the tangent vector of 
/, then paralleLtransport(/, v) represents a fixed direction with re- 
spect to that observer, so parallel transport is a physically important 
operation. 

Using just the functional definitions of the operations of geodesic 
tracing and parallel transport, it is easy to define more complicated 
functions which might be of interest, such as the following: 

paralleLcurve: {R ^ M) x V ^ [R ^ M), 
paralleLcurve(/, u) =7, 7: M ^ 7W, 

= geodesic(parallel_transport(/, (6) 

If we once again assume that the curve / is the geodesic world- 
line of a freely-falling non-rotating observer, and that v is orthogo- 
nal to the tangent vector of /, then, with respect to that observer, 
paralleLcurve(/, v) represents the world-line of an object which ap- 
pears stationary at a proper distance g{v, v) from the observer. [1] 

2.4.1 Computational accuracy and efficiency 

The benefits of 'lazy' function evaluation, mentioned in regard to vi- 
sualisation of functions in apply equally to numerical operations 
on functions. If, instead of being directly represented using functional 
programming techniques, functions are represented by their values at 
a finite number of points (as described in then a generic opera- 
tion on functions, such as that of parallel transport, would be forced 
to interpolate between (or, worse, extrapolate from) the known values 
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of the function / in order to obtain an estimate of f{t) for arbitrary 
values of t. Such a system is obviously both less accurate and less 
efficient than the evaluation of / only at the exact values of t required 
by the operation. 

3 Numerical experiments 

The analysis of the recent scientific claim j2| involved the modelling 
in GRworkbench of idealised interferometers orbiting the centre of our 
galaxy. The idealised interferometers each comprised a time-like curve 
representing the world-line of the beam-splitter of the interferometer, 
time-like curves representing the end-mirror of each arm of the inter- 
ferometer, and null geodesies connecting the beam-splitter world-line 
and the end-mirror world- lines, representing the world-lines of photons 
travelling along and back each arm of the interferometer. 

The world-line of the beam-splitter of each interferometer was 
modelled as a circular orbit in the equatorial plane of the spherically 
symmetric Schwarzschild black hole metric that was used to model our 
galaxy. Two different models for the world-line of the end-mirror of 
each arm were studied in the analysis of the recent claim. We do not 
describe them in detail here, but remark that in the more physically 
realistic of the two models, the world-lines of the end-mirrors are rep- 
resented by curves defined in a similar manner to the parallel curves 
of ^2A\ they lie at a fixed distance with respect to a non-rotating 
observer at the beam-splitter, although (unlike parallel curves) they 
do not lie in a fixed direction with respect to such an observer. (The 
interferometers modelled were in fact rotating very slowly.^) 

3.1 Geodesies defined by boundary conditions 

Given two nearby world-lines a and b, such as those of the beam- 
splitter and of one end-mirror, the null geodesies connecting them 
are physically important, as they represent the world-lines of light 
signals travelling between the observers represented by a and b. The 
problem of determining such null geodesies is an example of a 'two- 
point' boundary value problem, so called because there are boundary 

■^Each interferometer rotates at the same angular velocity as the beam-splitter's orbital 
motion about the centre of the gravitational field, which is the angular velocity of the 
Earth's orbital motion about the centre of our galaxy, which is very small. 
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conditions at two different values of the parameter of the null geodesic. 
For example, if we require that the intersection of the null geodesic 
/ with the world-line a be at a(0) (the 'photon emission' event), and 
that the intersection of the null geodesic with the world-line h be at 
/(I) (the 'photon reflection' event), then the boundary conditions on 
the null geodesic are specified as 



where tr is unknown until the boundary value problem is solved. Note 
that the requirement that the photon reflection event be at /(I) here 
(and in Equation below) is arbitrary; if /(I) = h{t^) then, by 
using a different affine parameter for the null geodesic /, we can make 
/(a) = 6(ir) for any other a > 0. 

There is another physically important definition of a geodesic in 
terms of a two-point boundary value problem: the unique geodesic / 
connecting any two nearby points in space-time. If the points are p 
and q then the boundary conditions are 



If / is time-like then it is possible for an observer to be present at both 
events represented by p and q; and if / is not space-like then the events 
represented by p and q are causally related, with p causally influencing 
g if / is future-directed, and the converse if / is past-directed. 

Because of the physical significance of these two types of geodesies 
defined by two-point boundary value problems, a general method for 
numerically determining them is an important tool for GRworkbench. 

3.1.1 Numerical determination 

To be concrete we describe here how GRworkbench determines geodesies 
defined by boundary conditions like Equation ((Hj). Later f ^:-{.1.2|) we 
briefly describe the similar procedure whereby geodesies defined by 
boundary conditions like Equation Q are determined. 
We require a function 



/(O) = a(0), /(I) = b{t,) 



(7) 



f{0)=p, f{l)=q. 



(8) 



connecting_geodesic : M y. M ^ {M. ^ M), 
connecting_geodesic(p, g) = /, /: M ^ 
/ is a geodesic satisfying Equation Q. 



(9) 
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GRworkbench implements connecting_geodesic by numerically deter- 
mining the tangent vector G Tp of / at p, and then returning 
geodesic (u) from Equation @. 

In determining the tangent vector v we are determining which 
direction, in space and time, to 'launch' a geodesic from p such that it 
'hits' q. This is accomplished in GRworkbench by minimising over u € 
Tp (using a standard multi-dimensional minimisation algorithm |5j) the 
amount A{u,q) by which the geodesic with tangent vector u at p 
'misses' the point q. 

We require of the 'miss' function A{u, q) that A{u, q) >0 and that 
A{u,q) =0 4^ geodesic(u) = connecting_geodesic(p, g), so that the 
desired global minima are simply local minima u such that A(u, q) = 0. 
Various definitions for A(u,q) are possible. (The particular choice 
of A{u,q) determines the 'basin of convergence', that is, the set of 
values of the initial guess x for which the numerical minimisation 
algorithm below, will converge to a local minimum such that 
^{u,q) = 0.) A simple definition is 

A{u,q) = ||0c(geodesic(u)(l)) -0c(g)||, (10) 

for some chart C, where ||-|| denotes the Euclidean norm on R". That 
is, the amount by which the point given by geodesic(u)(l) 'misses' the 
point q is defined as the Euclidean distance between the images of the 
two points on the chart C. 

The generic multi-dimensional minimisation algorithm, which at- 
tempts to numerically determine minima of a given function H : M" — > 
M, may be thought of as a function 

minimise: (M" ^ M) x ^ M", 

minimise(-ff, n*) = (a local minimum of H near u). (11) 

A minimisation over n S Tp is made possible by parameterising u by 
its components E M" on some chart C; the actual minimisation is 
performed over G M". 

Using the algorithm (|11|1 and the definition (|1()|) . GRworkbench 
is able to numerically determine the unique geodesic connecting two 
nearby points. A complete functional definition of connecting_geodesic 
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(continued from Equation Q) is 




— , v' = minimise(i7, (pdq) - (t^cip)), 



(12) 



where {4>c{q) ~ 'Pc{p)y ^ the coordinate direction from p to g on the 
chart C, has been used as the initial 'guess' direction in the numer- 
ical minimisation algorithm (|llj) . Like the 'miss' function described 
above, other definitions of the initial guess are possible, and the con- 
vergence (or otherwise) of the numerical minimisation operation will 
depend on the particular initial guess. If the space-time is sufficiently 
complicated between the two points then, like any numerical algo- 
rithm, the minimisation (|11|) may not converge to a solution u such 
that A(m, g) = 0, in which case the value of connecting_geodesic(p, g) 
is undefined. 

3.1.2 Determination of connecting null geodesies 

A slight modification of the algorithm described above is used by GR- 
workbench to determine null geodesies satisfying Equation 0. There 
are two important differences: (1) the vector u £ Tp must be null, 
so that the minimisation is to be performed over the null sub-space 
of Tp] and (2) the definition of the 'miss' function Anuii(w, ^) must 
refiect the requirement that the solution geodesic intersect the curve 
b of Equation ^ rather than the point q of Equation Q . 

The first difference is reconciled by minimising over three (non 
time-like) components (on some chart C) of the vector u, with the 
remaining (non space-like) component being determined by the re- 
quirement that u be null. The second difference is reconciled by using 
the following definition for the 'miss' function: 



where we have used A(u,q) from Equation pOj) . That is, the amount 
by which the point given by geodesic(M)(l) 'misses' the curve b is 
defined as the closest the curve gets to the point in the coordinates 
of a given chart. The one-dimensional minimisation over t in Equa- 
tion (|13|) is performed using a standard algorithm, [Sj implemented in 



/\nun{u,b) 



min A(u, 



(13) 
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Figure 2: An idealised interferometer simulated in GRworkbench, with 5 
orthogonal arms. The interferometer is orbiting the field centre, marked by 
the ball. The world-lines of the end-mirrors of each arm are joined to the 
world-line of the beam-splitter by null geodesies. 

GRworkbench with a functional interface like that of the minimisation 
algorithm with n = 1. 

4 Conclusion 

By 'gluing' together the methods described in this article, using their 
functional implementations, it is easy to define potentially complex 
numerical experiments. Figure |21 shows the interferometer described 
in Sjni simulated in GRworkbench. Physical properties of these inter- 
ferometer simulations, such as the travel time of photons as measured 
by the elapsed proper time along the world-line of the beam-splitter, 
yielded important results in our recent analysis. [21 As physical situ- 
ations continue to motivate the addition of new features in GRwork- 
bench, it becomes progressively more useful as a tool for numerical 
investigations in General Relativity. 



13 



References 



[1] Searle, A. C, GRworkbench, B.Sc. (Honours) thesis, Department 
of Physics, The Austrahan National University, 1999. 
Scott, S.M., Evans, B. J.K., and Searle, A. C, 'GRworkbench: A 
computational system based on differential geometry'. Proceedings 
of the Ninth Marcel Grossmann Meeting on General Relativity, ed 
Gurzadyan, V. G., Jantzen, R. T., and Ruffini, R., World Scientific, 
(2002) 458-467. 

Evans, B.J.K., Scott, S.M., and Searle, A. C., 'Smart Geodesic 
Tracing in GRworkbench', General Relativity and Gravitation, 34 
(2002) 1675-1684. 

[2] Moylan, A., Scott, S.M., and Searle, A.C., 'Can the Milky Way 

be weighed using Earth-based interferometry?', in preparation. 
Moylan, A., Numerical experimentation in GRworkbench, B.Sc. 
(Honours) thesis. Department of Physics, The Australian National 
University, 2003. 

For the original claim see Karim, M., Tartaglia, A., and Bokhari, 

A. H., 'Weighing the Milky Way', Glassical and Quantum Gravity, 
20 (2003) 2815-2825. 

[3] Moylan, A., Scott, S.M., and Searle, A. C., 'Functional program- 
ming framework for GRworkbench', in preparation. 

[4] For non-rotating observers who are not freely-falling, Fermi- 
Walker transport generalises the notion of a fixed spatial direc- 
tion; see for example Stcphani, H., General Relativity, Cambridge 
University Press, second edition, 1990, pp. 70-73. 

[5] Press, W. H., Teukolsky, S.A., Vetterling, W. T., and Flannery, 

B. P., Numerical Recipes in C, Cambridge University Press, second 
edition, 1992, pp. 412-418. 

[6] Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, 
B. P., Numerical Recipes in C, Cambridge University Press, second 
edition, 1992, pp. 402-405. 



14 



